THE EXPLICIT SIMPLIFIED INTERFACE METHOD 
FOR COMPRESSIBLE MULTICOMPONENT FLOWS. 



BRUNO LOMBARD* AND ROSA DONAT+ 

Abstract. This paper concerns the numerical approximation of the Euler equations for mul- 
ticomponent flows. A numerical method is proposed to reduce spurious oscillations that classically 
occur around material interfaces. It is based on the "Explicit Simplified Interface Method" (ESIM), 
previously developed in the linear case of acoustics with stationary interfaces (2001, J. Comput. 
Phys. 168, pp. 227-248). This technique amounts to a higher order extension of the "Ghost Fluid 
Method" introduced in Euler multicomponent flows (1999, J. Comput. Phys. 152, pp. 457-492). 
The ESIM is coupled to sophisticated shock-capturing schemes for time-marching, and to level-sets 
for tracking material interfaces. Jump conditions satisfied by the exact solution and by its spatial 
derivative are incorporated in numerical schemes, ensuring a subcell resolution of material interfaces 
inside the meshing. Numerical experiments show the efficiency of the method for rich-structured 
flows. 
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1. Introduction. Let us consider multicomponent flows composed of pure invis- 
cid fluids separated by material interfaces. These flows arise in a wide range of physical 
situations, from water-steam to bubbly flows, liquid suspensions or even high-speed 
impacts on solids. They may be modeled by the Euler equations, augmented by ad- 
ditional equations describing the fluid composition. The numerical simulation of such 
configurations leads to major difficulties. 

Indeed, even state-of-the-art numerical schemes for single-component flows cannot 
be applied directly in multicomponent flow simulations. These schemes give rise 
to oscillations and other computational inaccuracies near material interfaces. For 
example, any Godunov-type shock-capturing scheme which conserves the mass of 
individual species fails to maintain pressure equilibrium at a material interface yQ. The 
unphysical pressure oscillations stem from the fact that in a shock-capturing scheme, 
the transition across a material interface is governed by the numerical viscosity of the 
scheme. When two fluids are involved, intermediate states generated in the numerical 
transition layer (corresponding to the material interface) are not physically consistent 
with any component of the mixture: updating the pressure field via one particular 
equation of state generates erroneous pressure fluctuations. Since material interfaces 
lack the compressive mechanisms associated to shocks, errors generated within the 
diffused interface escape and contaminate all flow variables. 

Many numerical methods have been proposed to avoid these unphysical oscil- 
lations. See e.g. for a concise survey of up-to-date multicomponent methods. 
Remaining within the front- capturing (as opposed to front-tracking) computational 
framework, we distinguish two main approaches: 

Miscible models . Material interfaces are approximated by diffused fronts. An 
artificial equation of state (or mixture model) is defined, based on thermodynamical 
arguments [31 El EH > which is considered valid for the entire fluid mixture and is used 
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to update the pressure from conserved variables. These models have been analyzed 
in ^ |21 E3 Ht>l I31| . The elimination of unphysical oscillations generally involves 
sacrificing strict conservation to some degree, except in |'2,'i| where a conservative 
algorithm maintains oscillations under a computationally acceptable level. 

Purely inmiscible models, as in the present paper. A level-set function is used to 
track material interfaces in an Eulerian manner, i.e. without explicitely computing 
their location. Depending on the sign of the level-set, one applies the corresponding 
equation of state. Thermodynamic properties of the fluid change discontinuously 
across material interfaces, which are preserved as sharp discontinuities. A simple 
approach was first used in |24|. but it suffers from severe computational inaccuracies 
in multicomponent flow computations |16| . 

The Ghost Fluid Method (GFM) 8. is the best-known method belonging to this 
second group. On each side of a material interface, two fluids are considered: the 
"real fluid" (i.e. the fluid that really exists on this side) and a "ghost fluid" (i.e. 
the fluid with the same pressure and velocity than the real fluid on that side, but 
the entropy of the real fluid on the other side). Near a material interface, classical 
single-component schemes are then simply applied both on "real- fluid" values and on 
"ghost-fluid" values. Extensions to multidimensional problems and to other physical 
situations have been tackled by the GFM [7||^]. However, the GFM suffers from some 
inaccuracies which are not appreciated when material interfaces separate uniform 
states. These numerical artifacts are mainly due to the zeroth-order extrapolations 
used to define the ghost fluid. Simple minded first-order extrapolations do not work 

E 

The goal of the present paper is to remove these drawbacks of the GFM by cor- 
rectly increasing the precision of the computation across a material interface. We pro- 
pose to implement a carefully designed first-order extrapolation procedure to obtain 
the "ghost fluid" values. To do so, we adapt the Explicit Simplified Interface Method 
(ESIM), previously developed for the linear hyperbolic systems of acoustics with sta- 
tionary interfaces |21l 1221 12*7) . Note that another linear extrapolation in the context 
of the GFM has been proposed in 0] in order to couple Eulerian and Lagrangian 
computations. In the present paper, however, the goal is to enforce first-order jump 
conditions at material interfaces. 

Interface methods have been widely used in the numerical treatment of bound- 
aries and interfaces in PDE's: see e.g. the Immersed Interface Method (IIM), applied 
to elliptic equations ^H] an d linear hyperbolic systems [33]. See [2011221 f° r a concise 
survey of up-to-date interface methods. To our knowledge, the present paper is the 
first attempt to apply interface methods to the nonlinear hyperbolic system of the 
Euler equations. The computational complexity of the resulting algorithm is essen- 
tially the same as that of the GFM, and it can easily be coupled to a wide class of 
high-order shock-capturing schemes. 

The paper is organized as follows. Section 2 recalls the level-set framework to 
model compressible multicomponent flows. Section 3 describes the numerical schemes. 
The interface method is detailed in section 4. Numerical tests are proposed in section 
5. Conclusions and future works are drawn in section 6. 

2. The level-set framework for multicomponent flows. 

2.1. The Euler equations. We focus on multicomponent flows consisting of 
pure fluids separated by material interfaces. Assuming that all components can be 
described by a single velocity and pressure function, the flow can be modeled by the 
compressible Euler equations expressing conservation of mass, momentum, and energy 
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of the fluid mixture. Let p be the density of the fluid mixture, u the velocity, p the 



pressure and e = e + ^u 2 the specific total energy, with e the specific internal energy. 
Then, we have 



§-U + ^f(U)=0, (2.1) 
where the vector of conserved quantities U and the flux function / are 

f(U)= ( Pu 2 +P ] ■ (2.2) 
\ pe ) \ u(pe+p) J 

To close the system (|2.1f) . we need to specify the equation of state (EOS). In this 
paper, we consider the stiffened gas EOS 

P= (7- !)/>£- TPoo, (2.3) 

where p^ is the stiffness parameter. The basic polytropic gas case is recoverred for 
Poo = 0: in this case, 7 represents the ratio of specific heats. The EOS (|2.3|l is 
a reasonable approximation for gases, liquids, and even solids under huge pressure 
conditions 29 . The sound speed c is given by 
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Fig. 2.1. Two components Qo and f2i separated by a material interface. 

2.2. The level-set equation. The flow description is completed by an addi- 
tional equation describing the fluid composition. In level-set framework, a scalar 
function (j) is used to track material interfaces. For the sake of simplicity, we consider 
only two components f2o and fii, separated by a material interface at a(t) (figure 
I2.1|l . We suppose that each fluid satisfies l|2.3() , and that physical parameters of l|2.3(l 
may be discontinuous at a 

!(7o,Pooo) if x < a 
(2.5) 
(71, Pool) if x > a. 

The marker variable <p is initialized as the signed distance function to a(0), which is 
known initially. Hence its zero level-set defines the material interface at t = 0, while 
its sign determines the region occupied by each fluid. 

Since material interfaces propagate with the fluid velocity, the zero level-set of <j) 
identifies the material interface for all t > if 4>(x, f) satisfies the advection initial- 
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value problem 

86 86 

at dx (26) 

6{x,0) =x- a(0). 

The evolution equation in l|2.ti[) can be recast in conservation form, by combining it 
with the conservation of mass, 

Equations l]2.ip . I|2.3|l . and (|2.6|) (or $2.1\ ) form an "inmiscible" model in which the 
fluid "mixture" consists of either fluid £! or Qi, and where thermodynamic properties 
change discontinuously across the material interface. 

3. Numerical schemes. To integrate the Euler equations and the level-set 
equation for multicomponcnt flows, two strategies are available. The first one is 
to discretize the conservative system composed by the Euler equations (|2.1() together 
with the level-set equation in conservation form l|2.7[) . This discretization can be per- 
formed by applying the classical methodology of shock-capturing schemes for systems 
of conservation laws. It has been used in other papers (see e.g. [21]), but the addition 
of the level-set equation enlarges the Jacobian matrix, which increases the complexity 
of characteristic decompositions. This feature may become specially cumbersome in 
multidimensional calculations. 

The second strategy - simpler, and followed here - is considered in jjjj, where 
the Euler equations (|2.1|1 are discretized independently from the level-set equation 
in non-conservative form 12.6J1 . Following the Euler system is identified as the 
minimal system, hence it can be discretized with any standard conservative scheme, 
while the level-set equation is independently discretized following its advection form 
()2.6|l . The numerical flux function involved in the Euler equations will therefore not 
depend directly on 6, like in a single-component flow. The Lax-Wendroff theorem 
can still be applied ^J, hence the obtained numerical solution converges to a weak 
solution of the full conservative system. 

The discrete set-up to solve the system (|2.1|) - (I2.6|) is as follows: we consider a 
lattice of points in the x plane Xi = i Ax, where A a; is a uniform spacing parameter. 
Numerical values C/™ and 6™ are respectively considered as approximations to U and 
6 at Xi and time t n |3()| . To decouple the spatial and time integrations, we follow a 
method of lines approach ^1 f° r both the Euler system and the level-set equation. 
The semi-discrete approximation of the nonlinear hyperbolic system (|2.1|l is written 

4iUi=L a (U,i), (3-1) 
at 

whereas the semi-discrete approximation of the advection equation (|2.6|) is 

j- t 6 t = G{6,i). (3.2) 

The discrete operator Lq in (f 3 . 1 1> is specified in subsection 13. II We say that a grid 
point is regular if Lq uses numerical values belonging only to one fluid component. 
Otherwise, a grid point is called irregular, and the expression of Lq is modified by the 
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interface method detailed in section 4. The discrete operator G in i|3.2|) is specified 
in subsection 13. 21 

We emphasize the fact that the numerical schemes proposed in the next two 
subsections are by no means fundamental for the interface method. The reader's 
favorite single-component solvers can be adapted easily to the forthcoming discussion. 

3.1. Discretization of the Euler equations. When Xi is a regular point, all 
the values involved in the flux computations belong to one of the fluid components, say 
SI;. Then, the spatial discrete operator Lq = Lq 1 in i|3.1fl is written in the customary 
conservation form 

L Ql (U, i) = —£- (Fn, (Ui- s+ i, U i+S ) - F a , U i+S _i)) . (3.3) 

A x 

The numerical flux function F q in (|3.3|l is the trademark of the scheme; it is defined 
by a reconstruction procedure, and by a solver. In numerical experiments, we use 
ENO or WENO reconstructions E231E31. The width of the stencil s in J3~3f) is related 
with the theoretical order of accuracy of the spatial reconstruction (s = 3 for ENO-3 
or WENO-5). 

As a solver, we consider a flux-splitting construction |SJ^]. This choice has 
been considered in |23) within the mass fraction model for two-ideal gas flows. This 
flux-splitting requires two spectral decompositions of the Jacobian at each cell inter- 
face, which serve to perform upwind reconstructions of characteristic variables and 
fluxes. The additional cost is counterbalanced by the robust behavior of the scheme in 
pathology-prone situations. In addition, and as opposed to a Roe-type numerical flux 
function, no average-state needs to be computed at a cell interface, which is particu- 
larly useful for real-gas simulations. We only need to know the spectral decomposition 
of the Jacobian matrix: for the EOS l|2.3|l . this is given e.g. in |32| . 

The time integration of i|3.1|) is performed by the standard third-order TVD 
Runge-Kutta [HO], even when the fifth-order WENO-5 reconstruction is used. At 
follows from A x and from the classical CFL condition of stability 

CFL = . max (|<| + cf) ^ < 1, (3.4) 
z=0,..., N x Ax 

where N x + 1 is the number of grid points. 

3.2. Discretization of the level-set equation. We follow the same method 
of lines as for the Euler system. The spatial discretization is carried out using the 
Hamilton-Jacobi framework for the numerical approximation of the derivative terms 
|14l I26j . The upwind direction at Xi is determined by the sign of Uj. Full details are 
in Appendix A.l of 8 . 

Theoretically, solving (|2.6I) is sufficient to track the material interface. However, 
some extra-care must be taken for numerical purposes. Since u is nonuniform along 
the flow, the numerical representation of 4>{x, t) may become distorted |28| . leading to 
a poor estimation of the position of a. To keep <fi approximately equal to the distance 
function near the material interface, we follow a classical procedure proposed in |33j . 
and called reinitialization of the level-set. This procedure can be carried out in a 
number of different ways, here we follow [5| and solve to steady state the Hamilton- 
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Jacobi equation 



90 
~dt 



90 
dx 



-1=0, 



(3.5) 



t>(x,0) = (j)(x,t n ), 
where S is a smeared sign function given by 



S(0) 



The reinitialization equation is solved in fictitious time after each fully complete time 
step for the Euler equations, with a method of lines approach. For the spatial inte- 
gration of H3.5f) . we use a modification of Godunov's method detailed in Appendix 
A. 3 of 8 . For the time integration, we use the same TVD Runge-Kutta scheme as 
for the Euler equations. As a time step, we take At = Ai. After five integrations, 
we obtain a steady distance function, 0, which is then exchanged with 0. 



4. The ESIM for the Euler equations. 



4.1. The numerical interface treatment. The ESIM is a numerical treatment 
to be applied at irregular points, i.e. at points for which the flux computations in 
the right-hand side of H3.lt involve more that one fluid component. This treatment 
is carried out at each time integration (such as a Runge-Kutta substep). Let us fix 
an instant in time (to simplify the notations, the time variable is omitted from now). 
Then, the basic strategy of the ESIM can be schematically described as follows. 
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Fig. 4.1. Numerical values Ui and modified values U* . 



Consider the material interface at x = a separating the two media f^o and f2i 
C figure I2"T|) . We define smooth extensions U*(x) of the exact solution U(x): here, we 
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consider two-terms Taylor-like expansions past the interface, written 

d 

for x > a, U*(x) = U(a~) + (x - a) — U(a~), 

ox 

(4.1) 

fon<a, U*(x) = U(a + ) + (x-a)-^-U(a + ). 

ox 

Then, the interface treatment is divided in two parts. 

First, numerical estimations U* ofU*(xi) in (|4.1|) . called modified values (or ghost 
values, in the GFM framework), are sought at grid points surrounding a (figure l4~T)l . 
To do so, one must estimate U(a ± ) and -J^ U{ar c ). These last estimations, based on 
jump conditions presented in subsection 14.21 are detailed in subsection 14.31 

Second, at each irregular point Xi belonging to a particular fluid component, 
say f2;, the spatial operator L^i{XJ ,i) is constructed by considering only thermody- 
namically similar state values, that is: numerical values for grid locations in f2; and 
modified values at grid points on the other fluid component. This subject is detailed 
in subsection 14.41 

4.2. Jump conditions at material interfaces. Across a, the primitive vari- 
ables satisfy the classical zeroth-order jump conditions 

[u] = 0, [p] = 0, (4.2) 

where, for any function f(x,t), 

[/]= lim f(x,t)- lim f(x,t). 

x — >a+ x — >ol~ 

Since u and p do not jump across the interface, then their material - or Lagrangian - 
derivatives do not either. This fact immediately leads to the less-classical first-order 
jump conditions 



1 dp 
p dx 



0, 



2 ^ U 

PC d~x 



0. (4.3) 



One can also find second-order (and higher) jump conditions. However and unlike the 
acoustics case |20] > these jump conditions (not presented here) involve nonlinear 
combinations of p, u and their successive spatial derivatives, making the computations 
much more intricate. 

For polytropic gases, pc 2 = 'yp; since [p] = 0, the last equation of 1)4. 3|) can be 
simplified in [7 ^] = 0. Note also that nothing is said about [p] and [§^]. 

4.3. Numerical estimation of one-sided quantities. In a level-set model, 
the location of the material interface is given by a sign change in the marker variable 
4>. Suppose that 

0j x 4> J+ i < 0. (4.4) 

As a consequence, the interface lies somewhere between xj and xj+i. The subcell 
location of a is specified by the parameter 



(4.5) 
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which is readily estimated to second-order by a simple linear interpolation involving 
the smooth level-set function (f>. Indeed, defining 

« \4>j\ (46) 



\4>j\ + \<f> 



we have 8 = 6 + 0(Ax 2 ). Since the jump conditions Ij4.2(l and 14.3J1 concern the 
primitive variables W = T (p, u, p), it is easier to first look for estimations of W(a^) 
and ^ W{a^), and then to go back to the conserved variables U(a ± ) and ^ U(a ± ). 
In the following discussion and for any function /, the numerical estimation of f{a il ) 
or ^ /(a*) is denoted by / ± or /±. 

Density. The density is not subject to any constraint at a material interface, 
therefore numerical estimations p ± and p^ are performed via one-sided interpolations. 
On the left of a, elementary interpolations lead to 



(4.7) 



(4.8) 



p =-Up.,- 1 + {l + U)p.,, 

Px = (PJ - PJ-l) , 
while on the right of a, we get 

p+ = {2-e)p J+1 -(i-6) p J+2 , 

pi = (pj+2 - pj+x) ■ 

If p(x) is a piecewise C 1 function, we obviously have 

p(a ± )=p ± +0(Ax 2 ), ^{a k )=p£ + 0(Ax). 

Remark 1. The density estimations obtained above must remain within the 
physical range, i.e. p ± > 0. In some cases, negative values can be numerically 
obtained (e.g. in regions close to vacuum). Then, a different approximation procedure 
with constraints would be necessary (which is not considered here). 

Pressure. The pressure and its spatial derivative satisfy 1)4. 2|) and (|4.3|l 

p(a + ) = p(a~), 

i dp. i dp. ( 4 - 9 ) 



p{a + ) dx p(a ) dx 

From the definition of p and p^ , we have 

Pj = p~ + {xj-a)p~, 

(4.10) 

PJ+i = P +{xj+i-a)p+. 

We impose that numerical estimations satisfy exactly the same jump conditions as 
exact values. Combining 1)4. 11) [I with (|4.9|) leads to the 2x2 system 




1 (1 




(4.11) 
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The above matrix is always invertible, provided that p > 0. Inverting the above 
system gives the numerical values at a" 



P~ = ^+ ((l-6)^pj + dp J+1 ) , 

(1-0)^ + V p 
p 

Px = + \ (PJ+1 - Pj) 

(1-0) — + 



(4.12) 



To get estimations at a + , one uses <|4.12[1 and jump conditions (|4.9I) . Note that p + (or 
p~) is a convex combination with positive coefficents of p,j and pj+i, hence p ± > 0. 
As in the previous case, standard results allow to conclude that if both p(x) and p(x) 
are C 1 away from the material interface, then 

p(a ± )=p ± +0(Ax 2 ), ^(a ± )=p± + 0(Ax). 

Velocity. Zeroth-order and first-order jump conditions for the velocity are 

u(a + ) — u(a~), 

Applying the same procedure as for the pressure, we obtain a system which is invertible 
as long as {1 — 0)-f (p~~ +Pooo) /71 {p + + Poc 1) + > 0. This is always the case since 
p ± > 0. For stiffened gas EOS defining 

(4W) 

71 (p +P001) 

(for the ideal polytropic case £ = 70/71) with p~ given by (|4.12fl . we obtain 
U ~ = (1-0)^ + 6 (^-^^"J + ^^+i), 

(4.15) 

To get estimations at a + , one uses Ij4.15|l and the jump conditions (|4. lrifl . As before, 
under the appropriate smoothness assumptions, one has 

u(a ± )=u ± +0(Ax 2 ), ^(a ± ) = u I ± + 0(Ai). (4.16) 

ox 

Conserved variables. Once sided estimations and are computed, we 

deduce and from (|2.3|1 . These values are second-order approximations to the 
exact values under appropriate smoothness assumptions, hence 

!U + + (x - a) U+ for x < a 
(4.17) 
U + (x — a) U x for x > a 

satisfies U*(x) = U(x) + 0(Ax 2 ) near the interface. 
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4.4. Spatial integration at irregular points. From (|3.3fl and 14. 4|) . one de- 
duces that the irregular points are 

xjs+i, x J+s . (4-18) 

Modified values must also be computed at xj_ s and a;j+ s +i because the material 
interface can move between t n and t n+ \, crossing one grid point (but no more, be- 
cause of the CFL condition): as a consequence, these two regular points can become 
irregular points during one Runge-Kutta substep. The modified values are computed 
by substituting the appropriate grid point in the expresion (14.171) . Using l|4.(jf> as a 
second-order approximation to (|4.5|l . we obtain 

i = J+ 1,..., J+s + 1, U* =U~ + (i- J -6) AxU~, 

(4.19) 

i = J-s,...,J, U* =U + + (i- J -9) AxU+. 

At each irregular point Xi, Lq is now applied on numerical values on the same side 
than Xi, and on modified values on the other side than Xi, hence 

i = J — s + 1, J, 



L, t iU.r } = --^(Fao(U^ a+u ...,Uj,U* J+1 ,...,U* i+s ) 



(4.20) 



-Fno (Ui- S , Uj, U* J+1 , U*^^)) , 
i = J, J + s, 

Ln(U,i) = —±- (F ni (U*_ s+1 , U*j,U J+1 , U l+S ) 

-Fqi (U*_ s , U*j, Uj + i, U i+s -i)) . 

4.5. Summary of the implementation. Suppose that Runge-Kutta integra- 
tions (denoted by RK) have been performed up to the (m — l)-th substep. Then, the 
time-marching for both the Euler equations and the level-set equation can be summed 
up in 8 steps. 



Step 1 
Step 2 
Step 3 



location of the interface. Compute J (|4.4|) and 8 14.511 . 
construction of modified values. Compute l|4.19|l . 
construction of temporary values. Build two sets Ai and Bi 



{ C/f- 1) ,z = 0,...,J, f U*,i = J-s,...,J, 

A (m-i) = J B (m_1) = I 

\ U* ) i = J + l,...,J + s + l, 1 \ C/j m-1) , i= J + 1,...,N X . 

(4.21) 

Step 4: update Ai and Bi. Compute one Runge-Kutta substep of I3.1fl 
i = s,...,J+l, A\ m} = RK (a^jLqq^ , k<m 

(4.22) 

i = J, N x - s, Bj m) = RK (Bf \ L ni ) , k < m. 



Step 5: update ip. Compute one Runge-Kutta substep of l|2.6|l via u 



(m-l) 
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Step 6: location of the interface. Compute J l|4.4|) from 
Step 7: update U. From l|4.22|l and J, select 



m) 



f A[ m) for i = 0, J, 

{ B[ m) ioTi = J+l,...,N x . 



(4.23) 



Step 8: reinitialization of </>. If m = 3 (i.e. t — t n+ i), integrate (|3.5|l . 

4.6. Some remarks. Complexity. The algorithm presented in section 14.51 is 
simple. The computation of modified values U* does not depend on the discrete 
spatial operator. No analytical results - such as the solution of a Riemann problem - 
are required. As deduced from l|4.22l) . one does not need to write a new solver: the 
adaptation of known single-component solvers to the multicomponent case is direct. 
The only difficulty is to switch precisely modified values in the appropriate solver, 
depending on the sign of the level-set function, which is in fact an easy task. 

Computational cost. In comparison with single-component simulations, our ap- 
proach leads to a +25 % additional cost, both on a memory and computational time 
point of view. This cost is almost completely due to the level-set function (j>. Such a 
cost is inherent to level-set formulations, in order to know the composition of the fluid 
at each grid point. Eulerian methods that do not use level-sets require in counterpart 
to modify the Euler eigenstructure, leading to a similar additional cost. 

Consistency. As for the GFM the ESIM works because values used for time- 
stepping in (|4.22|l are thermodynamically similar: A* (respectively B*) in (|4.21|) 
satisfy the same equation of state. The algorithm amounts to consider separately two 
single-component flows, where no oscillations exist. 

Single- component flow. In the limit case 70 =71, Pooo = Pool, the flow is single- 
component, and material interfaces amount to classical contact-discontinuities. The 
ESIM behaves equally well, sharpening contact discontinuities. 

Conservativity. The ESIM is formally non-conservative locally, since 



like in or [H]. However, the lack of conservation is only introduced at one cell 
boundary on the entire domain. In practice, this feature does not seem to spoil 
convergence to the correct solution. 

Note that a fully-conservative version of the GFM has been developed [35] . Such 
an approach could probably be adapted to the ESIM, but we do not look further in the 
present paper: unlike in |25| where inert shocks and detonation waves are adressed, 
we only focus on the material interfaces (where conservation errors are not crucial). 

5. Numerical experiments. 

5.1. Configurations. Four numerical experiments are proposed. Tests 1 and 2 
illustrate the interaction of shock waves with one and two material interfaces. Test 
3 is a pure advection problem, with smooth structures. Test 4 concerns nonlinear 
acoustics. Analytical values and numerical values are respectively shown in solid lines 
and dotted lines, and we take CFL=0.66. 

All tests have initially-isolated material interfaces. Riemann problems are not our 
concern, since jump conditions l|4.2|l - which are a building-block of the ESIM - are 
then generally not satisfied. Note that the GFM has the same limitation. In practice, 



Fn (Uj. s+1 , Uj, U* J+U U* J+S ) ? F na {U \ 



J-S+l! 



— ) Uj, C/j+i, U j +s ) , 

(4.24) 
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the GFM often works for Riemann problems, but some failures (see e.g. data of Test 
4 in [2]) can be explained by the fact that basic assumptions are not satisfied. 



Density Velocity 




0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 I 0.1 0.2 03 0.4 0.5 0.6 0.7 0.8 0.9 



Fig. 5.1. Test 1-a: Mach 1.95 shock wave interacting with a material interface. WENO-5 
coupled with the ESIM (exact values: solid line; numerical values: points). 

5.2. Test 1: shock-interface interaction. Firstly, we consider a shock-interface 
interaction problem, previously studied in |ril| . On aim long domain, the initial con- 
figuration consists of a stationary material interface at ao — 0.5 m, separating two 
fluids with different EOS: a polytropic gas on the left, a stiffened gas on the right. A 
left-going Mach 1.95 shock wave is initially set at ot\ = 0.7 m. Physical parameters 
are 

po = 1.000 kg/m 3 , po — 1 Pa , u = m/s ,70 = 1.4, poo = Pa, 
< pi = 5.000 kg/m 3 , pi = 1 Pa , u± = m/s ,71 =4, p x 1 = 1 Pa, 
. p 2 = 7.093 kg/m 3 , p 2 = 10 Pa ,u 2 = -0.7288 m/s , 72 = 4.0, Poo2 = 1 Pa. 

(5.1) 

At t = 4.69 10 2 s, the shock wave collides with the material interface, leading to a 
left-going shock, a left-going material interface, and a right-going rarefaction fan (see 
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figure 4 in [3T]). 

Density Velocity 




0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0.1 0.2 0.3 0.4 0.5 0.6 0.7 OS 0.9 I 

Km} x(m| 



Fig. 5.2. Test 1-b: Mach 5 shock wave interacting with a material interface. WENO-5 coupled 
with the ESIM (exact values: solid line; numerical values: points). 

Figure [3TTI shows exact values and numerical values of p, u, p, and entropy S at 
t = 0.202 s (after 200 time steps). The computations are performed on N x — 200 grid 
points with WENO-5 coupled with the ESIM. The agreement between analytical and 
numerical values is good. The main interest of this example is to show that the ESIM 
is robust and behaves well, even when a shock wave is in the vicinity of the material 
interface. No spurious oscillations induced by Taylor expansions are seen in u and p 
around the material interface (near x = 0.32 m). 

Secondly, we consider the case of a stronger shock wave. The physical parameters 
around the material interface are the same than in Q5.1[l : only the parameters in the 
post-shock fluid are modified 

p 2 = 8.116 kg/m 3 , p 2 = 77.80 Pa , u 2 = -2.428 m/s , 72 = 1.4, Poo 2 = 1 Pa, 

what amounts to a Mach 5 left-going shock wave. The wave phenomena (after the 
collision between the shock wave and the material interface) arc the same than in the 
previous example. The computations are performed with WENO-5 coupled with the 
ESIM. 
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Figure I5"2*l shows exact values and numerical values of p, it, p, and entropy S at 
t = 0.112 s (after 300 time steps). The agreement between numerical and analytical 
values is good, even for this strong shock test. No oscillations are visible in u and p 
around the material interface (near x = 0.26 m). 



Density 



Velocity 



I 



0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 




0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 



Entropy 



Pressure 



0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 



0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 



Fig. 5.3. Test 2: shock-bubble interaction. WENO-5 coupled with the ESIM (fine grid solution: 
solid line; numerical values: points). 

5.3. Test 2: shock-bubble interaction. We consider the ID version of the 
classical shock-bubble interaction problem, numerically addressed e.g. in [25 . with 
a conservative algorithm. On a 0.445 m long domain, a helium domain, initially 
delimitated by a = 0.2 m and «i = 0.25 m, is at rest in air. A left-going Mach 1.22 
shock wave is initially set at a-i — 0.275 m. Physical parameters are 



Pi) 


= 1225 kg/m 3 , po = 


1.01325 10 5 Pa , M = 


0.0 m/s ,70 = 


1-4, PooO = 


Pa, 


pi 


= 0.2228 kg/m 3 , Pl 


= 1.01325 10 5 Pa , iti 


= 0.0 m/s , 71 


= 1.648, poo 


i = Pa 


P2 


= 1225 kg/m 3 , p 2 = 


1.01325 10 5 Pa , u 2 = 


0.0 m/s ,72 = 


1.4, Poo2 = 


Pa, 


P3 


= 1686 kg/m 3 , P3 = 


1.59059 10 5 Pa , u 3 = 


—3.59 m/s , 73 


= 1.4, Poo3 


= Pa. 
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See [23 for a description of wave phenomena. Since no analytic solution is available, 
we compute the solution with WENO-5 coupled with ESIM on a fine grid with 3200 
mesh points; we refer to this approximation as the "exact solution" for comparison 
purposes. 

Figure shows "exact values" and numerical values of p, u, p, and S at t = 
2.86 10~ 3 s (after 400 time steps). The computations are performed on N x — 400 grid 
points with WENO-5 coupled to the ESIM. The agreement between "exact values" 
and numerical values is very good. The resolution is better than in |231 : compare the 
density inside the helium in figure IB~3*I with the density in figure 2 of 23 . 

To conclude the tests 1 and 2, let us notice that similar results could be obtained 
with the GFM treatment (not shown here). Indeed, when only flat profiles are in- 
volved, the GFM and the ESIM have a similar behavior. When a shock wave collides 
with a material interface, both treatments produce a "sloping" behavior in density 
and entropy very close to the "overheating" phenomenon found in shock reflection 
problems |H| ■ Increasing the Mach number seems to accentuate this " overheating" , 
producing in addition slightly perturbed post-shock values: small acoustic perturba- 
tions can be observed downstream in figure l5~2l More testing (not shown) has been 
performed with Mach numbers up to 9, which confirm these observations. No crash 
due to negative pressure values has been observed on these tests. 

The goal of the next two tests is to show the advantage of using the ESIM (com- 
pared to the GFM) for rich-structured flows. 



Method 


N x 


L\ error 


L\ order 


Method 




L\ error 


L\ order 




100 


6.44e-3 






100 


3.76e-3 




ENO-3 


200 


1.62e-3 


1.99 


ENO-3 


200 


4.74e-4 


2.98 


+ 


400 


6.86e-4 


1.24 


+ 


400 


6.42e-5 


2.88 


GFM 


800 


2.40e-4 


1.51 


ESIM 


800 


8.62e-6 


2.89 




1600 


6.44e-5 


1.89 




1600 


9.50e-7 


3.18 




3200 


2.05e-5 


1.65 




3200 


1.21e-7 


2.98 



Table 5.1 
Measures of convergence in Test 3. 



5.4. Test 3: pure advection. We consider a 1 m long domain with two material 
interfaces initially at a® — 0.160 m and a± = 0.526 m. The pressure and the velocity 
are initially constant: u(x, 0) — 300 m/s, p(x,0) = 10 5 Pa. The density is initially 



p(x,0) 



1 + 0.3 sin(50 (x — «o)) kg/m 3 if ao < x < a%, 



1 kg/m else. 
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Density (GFM) Density (ESIM) 




Entropy (GFM) 




Pressure (GFM) 




Entropy (ESIM) 




0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 



Pressure (ESIM) 



0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 



0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 



Velocity (GFM) 



Velocity (ESIM) 





Fig. 5.4. Test 3-a: 200 grid points. ENO-3 coupled with the GFM (left column) and with the 
ESIM (right column). Exact values: solid line; numerical values: points (note that the scales for p 
and u are magnified). 
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Density (GFM) 



Density (ESIM) 




0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.S 0.9 




0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 



Entropy (GFM) 




Entropy (ESIM) 



It mm 



0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 




0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.S 0.9 



Pressure (GFM) 



HKIKKtp 
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100060 ■ 
100040 ■ 
100020 - 
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»(i - 
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Pressure (ESIM) 



[(M)HH) 
100080 
100060 
100040 

km]u:o 

lUflUN) 

99980 
99960 
99940 
99920 



0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 



0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.S 0.9 



Velocity (GFM) 



31)0.5 r- 
300.4 - 



299.9 ■ 

2W.K - 
299.7 - 
299.6 ■ 
299.5 L 



Velocity (ESIM) 



' 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.S 0.9 




Fig. 5.5. Test 3-b: 800 grid points. ENO-3 coupled with the GFM (left column) and with the 
ESIM (right column). Exact values: solid line; numerical values: points (note that the scales for p 
and u are magnified). 
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Method 


N x 


A(p) 


Order 


A(pu) 


Order 


A(pe) 


Order 




100 


7.25e-l 




2.21e+2 




3.22c+4 




ENO-3 


200 


3.82e-l 


0.92 


1.14e+2 


0.95 


1.72e+4 


0.90 


+ 


400 


2.21e-l 


0.79 


6.65e+l 


0.77 


9.97c+3 


0.78 


GFM 


800 


1.27e-l 


0.80 


3.83e+l 


0.79 


5.74e+3 


0.79 




1600 


7.10e-2 


0.84 


2.13e+l 


0.84 


3.19e+3 


0.84 




3200 


3.61e-2 


0.97 


l.Q8e+l 


0.98 


1.62e+3 


0.97 




100 


1.98e-l 




5.96e+l 




8.94c+3 




ENO-3 


200 


4.69e-2 


2.08 


1.40e+l 


2.09 


2.11c+3 


2.08 


+ 


400 


1.07e-2 


2.13 


3.23e+0 


2.11 


4.85e+2 


2.12 


ESIM 


800 


2.75e-3 


1.96 


8.27e-l 


1.96 


1.24c+2 


1.96 




1600 


6.82e-4 


2.01 


2.37e-l 


1.80 


3.41c+l 


1.86 




3200 


1.54e-4 


2.14 


6.72e-2 


1.81 


9.36c+0 


1.86 



Table 5.2 
Conservation errors in Test 3. 



Physical parameters are 



(7, Poo) 



70 = 1.40, pooo = 10 4 Pa if x < a , 

71 = 1.67, pod = 10 5 Pa if ao < x < a\, 



72 = 1.40, poo2 = 10 4 Pa if x > a x . 



This configuration amounts to an advcction equation for p; material interfaces are 
advected at the velocity u; p and u remain theoretically constant. Because of the 
discontinuous physical parameters, the entropy S is discontinuous at ao and at ot\. 
Numerical experiments are performed with ENO-3, coupled with the GFM or with 
the ESIM. 

Figures 15.41 and 15.51 show exact values and numerical values of p, S, p, and u 
at t = 6.62 10~ 4 s, respectively for N x — 200 grid points (hence 25 grid points by 
wavelength) and N x = 800 grid points (hence 100 grid points by wavelength). Notice 
that the zeroth-order extrapolation of S used by the GFM result in jumps of p at 
material interfaces, that are advected with the flow; these glitches are also transferred 
to u, and p, near x = 0.35 m and x — 0.75 m, and act as sources of acoustic noise. 
These glitches still exist with 800 grid points ( figure 153) 1 with the GFM. On the other 
hand, no entropy or density glitches are observed when the ESIM is used instead. 

We must mention that ENO reconstructions of a sinusoidal profile might produce 
oscillations on the level of the truncation error. These small spurious oscillations are 
seen in p and u, which should have flat profiles, and can be observed even without 
material interfaces. When these oscillations interact with the GFM-produced glitches 
at material interfaces, they are amplified and lead also to spurious acoustic waves. 
These oscillations are entirely reconstruction dependent and, since they are of the 
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order of the truncation error (see also our convergence measures) , they are essentially 
not visible on 800 grid points when the ESIM is used (figure IS3)l . 

Measures of convergence are provided in table 5.1. Coupled with the ESIM, 
the ENO-3 scheme maintains third-order convergence. Coupled with the GFM, the 
ENO-3 scheme looses accuracy and shows a 1.6 order of convergence. Table 5.2 shows 
measures of conservation errors induced by the GFM and by the ESIM, for various 
values of N x (see subsection 14 .6f) . To do so, we compute 



AU(T 1 N X ) = max 

n=0,...,N t 



^- u ") +n ^nut)-f{uD) 



(5.2) 



at T = 1.05 10" 3 s, where / is the flux function , and N t = Trunc {T/ At). 

The domain of measure is bounded by io — 10 and i\ = N x — 10. Getting fiable and 
meaningfull measures requires some care. Indeed, errors of conservativity vary a lot 
with the position of a inside the meshing, hence with t. The "max" in (|5.2|l ensures 
almost-steady values of A U(T, N x ) and reliable measures. 

A fully conservative scheme would satisfy A U(T, N x ) = for all values of T and 
N x (see e.g. (12.33) in [El)- Here, errors are non-null, but they decrease with Ax. 
Conservation errors induced by the ESIM are much smaller than that induced by 
the GFM. The table 5.2 indicates first-order conservation errors for the GFM, and 
second-order conservation errors for the ESIM. 

Measures of convergence and conservativity have been done also with the WENO- 
5 scheme coupled with the GFM or with the ESIM. The results are essentially the 
same than in tables 5.1 and 5.2. To get a fifth-order accurate scheme (as in [37]) and 
a higher order of conservation would require to develop a higher order ESIM, which 
is not investigated in the present study. 



1.0015 



10 20 30 40 50 00 70 30 00 100 

"(ml 

Fig. 5.6. Initial values of the pressure p at to = 5.1 10 -2 s in Test 4- 

5.5. Test 4: nonlinear acoustics. As a last example, we consider a rich- 
structured wave interacting with a stationary material interface. This example is 
a generalization of the Test 1 studied in |27J : for small amplitude initial data, the 
simple linear equations of acoustics used in |27j are valid, and the wave remains 
smooth. Higher amplitudes of initial data lead to a nonlinear problem, where shocks 
can develop. 
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Fig. 5.7. Test 4. s = 10~ 3 (left column), e 
of fine grid solution (solid line) and numerical 
Errors with the ESIM: (c), (d). Errors with the 



= 1 (right column) at t± = 8.6110 2 s. Snapshots 
values (dotted line) of p with the ESIM: (a), (b). 
GEM: (e), (f). 
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To illustrate both cases, we use exactly the same physical parameters than in 
|27| . On a 300 m long domain, a material interface separates piecewise constant 
values, initially motionless 

po = 1000 kg/m 3 , po = 10 5 Pa , u = 0.0 m/s , 7o = 3, p^o = 7.50 10 8 Pa, 

p x = 1200 kg/m 3 , p! = 10 5 Pa , u x = 0.0 m/s ,7i = 4, Poo i = 2.35 10 9 Pa. 

From l]2.4p. one deduces sound speeds: c = 1500 m/s, and ci = 2800 m/s. These phys- 
ical parameters correspond respectively to water and Plexiglass under atmospheric 
pressure. A perturbation AWo(i) is added to initial data on medium Oo, which is 
an exact solution of acoustics and leads to a right-going wave 

A W (x) = ef(t --) T (-% —52-, -po) . (5.3) 
V Co J V c o Po c o J 

The function / is a C 5 spatially-bounded sinusoid 



q 1 
^a fe sm(/3 k w c £) if0<£<— , 



■■■{ *=i /c (5.4) 

else, 

with /3fe = 2 fc ~ 1 , w c = 27r/ c ; the coefficients afc are: a\ = 1, 02 = —21/32, 03 = 63/768, 
04 = —1/512. The central frequency is f c = 50 Hz. Elementary calculations show 
that 

„ „ „ e . (5.5) 
Po c Po 

For e = 0, the solution remains a stationary material interface. For small values 
of e (typically e = 10~ 3 ), the acoustics limit is valid and one obtains a pure right- 
going wave. For higher values of e (typically e = 0.1), the perturbation (|5.3|) does 
not satisfy exactly the nonlinear Eulcr equations: the wave separates into a weak left- 
going wave and a right-going wave. In both cases, the right-going wave is reflected and 
transmitted by the material interface. The analytical solution is not detailed in the 
linear case |20j . No analytical solution is available in the nonlinear case: one computes 
an "exact solution" on a fine grid of 3200 grid points. Numerical experiments are 
performed with N x = 400 grid points, which is 40 points in the wavelength A c = Co// c , 
and to = 5.1 10~ 2 s. Initial values of the p at to are shown in figure 15^1 for e = 10 -3 . 

Figure ISTI shows results at t\ = 8.61 10 -2 s, after 200 time steps. Left and right 
columns concern respectively e = 10~ 3 and e — 0.1. Figures [5.71 fa) and (b) show 
numerical values and exact values of p computed by WENO-5 coupled to the ESIM. 
Logically, (a) is similar to figure 3-(d) of 27 . In the nonlinear case (b), the reflected 
and transmitted waves have developed a shock. In both cases, the agreement between 
exact values and numerical values is excellent. 

To see clearly differences between the ESIM and GFM treatments, we display 
errors for both methods: with the ESIM in (c) and (d), with the GFM in (e) and (f). 
In the nonlinear cases (d) and (f), errors are displayed from x = 85 m to x = 120 
m, to avoid the shock area. In both cases, the graphs confirm that the ESIM is more 
accurate than the GFM. 
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Unlike in |27], the material interface is allowed to move with the flow. For 
e = 1CP 3 , the measured movement is lower than 10 -7 m: so, the approximation 
of stationary material interfaces is justified. For e = 0.1, the material interface moves 
from 96.3 m to 96.114 m, and then it comes back to the initial position 96.3 m: it is 
logical, since the velocity is symetric with respect to zero. 

6. Conclusion. We have proposed an extension of the " Explicit Simplified Inter- 
face Method" (ESIM), previously developed in acoustics to treat material 
interfaces in ID multicomponent Euler flows. 

The method enforces the numerical solution to satisfy zero-order and first-order 
jump conditions at the material interface and can be coupled with the user's favourite 
high-order shock-capturing scheme for single component flows. It behaves as robustly 
as the GFM [5] in numerical simulations involving flat states, while it displays a 
superior performance in flows with rich structures. 

Our numerical simulations show that the interface treatment we propose guaran- 
tees a numerical solution with no unphysical numerical artifacts due to the material 
interface. The behavior observed for the numerical errors is similar to that observed 
for the same underlying method when applied to a flow without material interfaces. 

This paper was focused on ID flows: applying the ESIM to 2D cases is a challeng- 
ing project, subject of future works. Some key tools for the 2D implementation have 
been validated already in the simple case of 2D linear acoustics with stationary inter- 
faces 122- Other ingredients (such as multidimensional level-set-based extrapolations) 
may be found in [SUT^]. 
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